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. Abstract 

CJ ■ 
• 1— I . 

We derive analytic expressions for infinite products of random 2x2 matrices. The 
(-H , determinant of the target matrix is log-normally distributed, whereas the remainder is a 

^ l' surprisingly complicated function of a parameter characterizing the norm of the matrix 

and a parameter characterizing its skewness. The distribution may have importance 
I ' as an uncommitted prior in statistical image analysis. 

> ■ 

' 1. Introduction: Considerable effort has been invested over the last half century 

. in determining the spectral properties of ensembles of matrices with randomly chosen 

' elements and in discovering the remarkably broad applicability of these results to 

, systems of physical interest. In spite of a similarly rich set of potential applications 

' (e.g. in the statistical theory of Markov processes and in various chaotic dynamical 

, systems in classical physics), the properties of products of random matrices have 

' received considerably less attention. See ref. for a survey of products of random 

O , matrices in statistics and ref. ||| for a review of physics applications. 

c/3 ' The purpose of the present manuscript is to consider in some detail the limit for 

^ oo of the ensemble of matrices 

' where t > is a real parameter and the X„ are real d x d matrices with all elements 

^ , drawn at random on a distribution of zero mean and unit variance. If this distribution 

has compact support, the probability that the matrix Y should become non-positive 
definite vanishes for N —^ oo. In one dimension, d = 1, it is well-known from the 
law of large numbers that log Y has a Gaussian distribution, but because of the non- 
commutativity of matrix products, the distribution is much more complicated for 
d>2. 

In this paper we derive some general properties for the limiting distribution V{Y) 
and determine it explicitly for d — 2. In section |^ we establish a compact diffusion 
equation for the distribution valid for any d. In section ^ we derive a simple expression 
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for any average over the distribution, and we show that the determinant det [Y] has 
a log-normal distribution. Sections ^ and ^ will be devoted to the determination of 
the explicit form of V for d = 2. We shall first write the diffusion equation using 
an appropriate parameterization of Y . The resulting partial differential equation will 
then be solved subject to the boundary condition that ViY) supports only the identity 
matrix in the limit of r ^ 0. This explicit solution will require new integrals involving 
Jacobi functions. The derivation of these integrals will be given in the Appendix. 

2. The diffusion equation: The normalized probability distribution is (for 
given N and variable r) 



Vn{Y,t) 



Y 



N 

n 
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Xi,... ,Xn 



where the integrand is a product of (5-functions for each matrix element of Y and the 
average runs over all the random matrices. Pealing off the A^th factor in the product 
and using only that the A"„ are statistically independent, we derive the following exact 
recursion relation 



VNiY,T) = (dct 



N 



Yil 
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where the average is over the iVth random matrix, here renamed X. The determinantal 
prefactor of Vn-i is the Jacobi determinant arising from the general matrix rule 



5[Y - ZM] - 
with M = 1 + y/r/N X. Since 

d{ZM\ 



5[YM-^ - Z\ 



det 



d(ZM) 
dZ 
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dZ, 
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the Jacobian is block diagonal with d identical blocks, and the prefactor follows. 

The recursion relation (||) is of the Markovian type with the initial distribution 
'Po{Y,t) = S[Y — 1]. It converges for iV — > oo under very general conditions (which 
we shall not discuss here) towards a limiting distribution V{y, r) = limAr^oo VNiy, t). 
Expanding the recursion relation to O (l/N) and using the fact that all the matrix 
elements of X are statistically independent with zero mean and unit variance. 







we obtain to leading order 

+ id + 
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with implicit summation over all repeated indices. The assumed convergence towards 
a limiting distribution requires the expression in the parenthesis to vanish in the limit, 
so that 

This is a diffusion equation of the Fokker-Planck type with r playing the role of time. 
It must be solved subject to the initial condition that Viy, 0) = S[Y — 1]. 

Both the diffusion equation and the initial condition are invariant with respect to an 
orthogonal transformation Y AI^YM, where M is an orthogonal matrix satisfying 

M = 1. Since the number of free parameters in an orthogonal transformation is 
^d{d — 1), the number of "dynamic" variables in the distribution is d? — ^d{d — 1) = 
^d{d + 1). Since the distribution only has support for det [F] > 0, this number 
consists of d independent eigenvalues and ^d{d— 1) rotation angles in a singular value 
decomposition. 

For d = 1 the solution to (|^) which approaches (5[y — 1] for r ^ is 



rd=iiY)^—^cxp 



(logy + r/2)^ 
2t 



(8) 



As expected, it is a log-normal distribution. 

3. Averages: Remarkably, equation may be written in the much simper form 

(9) 



dV _ IdHY.kY.kV) 



dr 2 dYudYji 

without any explicit reference to d. Defining the average of a function f{Y) by 

(/> = j f{Y)V{Y)dY (10) 
with dY — dYij , we obtain from (^) 

'''' A(y>^.JJ^ ■ (") 



dr 2 \ ' dY,edYje , 

This equation permits in principle the determination of the moment of any product 
of matrix elements. The first two are found to be 

{Y,) = 6., (12) 
(Y^jYki) = e-'^S.kSji (13) 

The exponential growth of the averages with "time" r is a consequence of the multi- 
plicative nature of the problem. 

The determinant D = det [Y] is, according to the definition of the product (^, 
an infinite product of random real numbers that converge towards unity, and log I? 
must have a Gaussian distribution according to the law of large numbers. Its mean 
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and variance are, however, different from those of the one-dimensional distribution 
(P). The distribution of the determinant is also an average 



F{D) = {S{D-dct[Y])) . (14) 



Using the fact that 



(15) 



we obtain the following equation for F 



Apart from the factor d in front, this is identical to the diffusion equation (^) in one 
dimension. Consequently the determinant has a log-normal distribution 



F{D) = ^= exp 



{\ogD + Td/2f 

2Td 



(17) 



which is obtained from (|^) by replacing r by rd. The distribution has support only 
for positive values of D. It can be shown in general (and we shall demonstrate it 
explicitly for d = 2 below) that the distribution of the determinant factorizes in V. 

4. The case d — 2: The first non-trivial case is d = 2 where the general matrix 
is first parameterized using a quaternion or 4-vector notation 

In this representation the determinant becomes a metric with two "space" and two 
"time" dimensions 

D = Y,^^ Yl + Yi - Yi . (19) 

The structure of this expression and the positivity of D suggest the following param- 
eterization in terms of one imaginary and two real angles 

Yq = \/Z? cosh -0 cos e (20a) 

Yi = \/Z? sinh ip cos (j) (20b) 

Y2 = ^ cosh ip sin d (20c) 

/d sinh ?A sin . (20d) 

The Jacobi determinant of the transformation from {Yo, ^i, i^2, ils} to {D,ip,9,(j)} is 
simply 

J- DsinhV'coshV' . (21) 
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Orthogonal 2x2 matrices are generated by the matrix ( " ) , which is associated 
with Y2. Thus, an orthogonal transformation rotates the angle (j), and V{Y,t) must 
be independent of as indicated above. 

In these variables the diffusion equation simplifies to 



dr ' " dD ' 9D2 

I d^V 1 dV 1 ff^V 

+ -{1 + tanh^ + -(tanh V. + coth^)— + . (22) 

Taking into account the factor of D in the Jacobi determinant, we replace the original 
distribution V with the product of the determinant distribution F{D) given in (17) 
and an as yet unknown function of V' and 0, 

V = ^F{D)G{4>.e) , (23) 

and find that G satisfies the diffusion equation 

- . -(1 + tanh^^)— + -(tanh^ + coth^)- + -— . (24) 
The corresponding normalization integral is found from the Jacobi determinant, 

p2-K poo 

I d0 dtp 2 sinh 2^- G(V', 9) ^ I . (25) 
Jo Jq 

This normalization integrals ( p5| ) suggests that it is more convenient to employ still 
another variable 



+ Y,^ + Yj + 
Yi-Y^+Yi-Y^ 



cosh 2^ = "2 '2 . J2 J2 ■ (26) 



With this variable the normalization integral takes the form 



2tt 



de / dzG{z,e) = 1 , (27) 



and the diffusion equation ( p^ ) becomes 

dG 1 2z d^G dG , . ,d^G 

dr Az + l dO^ dz ^ ' dz^ ^ ' 

This equation must be solved with the boundary condition that V{Y) in the limit 
T reduces to a product of delta functions which select only the identity matrix. 
This evidently requires Yq 1 and Yi,2,3 and, consequently, — > 1, z ^ 1, and 
— > 0. Since F{D) — > 6{D — 1), the initial condition takes the form 

G{z,0)^S{z-l)S{9) (r^O). (29) 

The limiting distribution should be approached from above (i.e. from z > 1). 
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The form of the diffusion equation 1^8^) reveals that G may naturaUy be expanded 
in a Fourier series 



oo 

^(^'^) = ^ E G„(z)e™« (30) 

with coefficients that obey 



2tt 

n— — OG 



^ - "i" TTT^" + ' ^'^^ ■ ^^^^ 

For the special case n = 0, we recognize Legendre's differential operator on the right. 
The normalization condition only affects Go and becomes 

dzG'o(z) = l. (32) 

The initial condition (|29|) implies that 

Gn{z) Siz ~ I) (r^O) (33) 

for all n. 

5. Explicit solution: All that remains is to determine the angular functions 
Gn{z). One relatively simple way is to use Sturm-Liouvillc theory, and we now outline 
the main steps in this procedure. 

The differential operator ( "Hamiltonian" ) appearing on the right hand side of 
eqn. (|3l| ) may be written 

^ , d , ^ . d n'^ 2z 

^=^^ -Ih 7' (34 

dz^ ' dz A z + 1 ^ ' 

which shows that it is Hermitean. Let the spectral variable (which dcnumerates the 
eigenvalues and may be both discrete and continuous) be denoted r, and let gn\z) 

(r) 

be the eigenfunction corresponding to the eigenvalue A„ , 

ngi^Hz)^\i:^gi:\z) . (35) 

The Hermiticity of Ti. guarantees that the eigenvalues are real and that the eigenfunc- 
tions arc both orthogonal and complete on the interval 1 < z < oo, 

/ dzg(:Hz)gl:'\z)^'-^ (36) 

Y,fi^:^gi^\z)gi:Hz')^Siz-z') (37) 

r 

(r) 

with a suitable measure, /i„ . 

The solution of the diffusion equation (|T]) with initial condition ( p3| ) takes the 
form 

G„(z, r) = J2 ^^^r[^g^'Hl)g^nHz) cxp (A(:)r) . (38) 
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In view of the completeness (p7|), these functions indeed satisfy the initial conditions 
at r = 0. The appearance of ^^'^''(l) in this expression requires the eigenfunctions to 
be regular at z = 1. 

We now present the complete solution of the eigenvalue problem. (Further details 
are given in the Appendix.) The eigenvalue spectrum contains discrete values (for 
n> 2) as well as a continuum 



1^2 



\-\n^-\~t^ 0<t<oo. 
The properly normalized discrete eigenfunctions are Jacobi polynomials 



n/2 



' (40) 

while the eigenfunctions in the continuum are Jacobi functions of complex index 

/ , \ n/2 

- (4^) p^^:U;.,j^) (41) 

with the measure obtained from the integral ( |3^ ) as 

\tta,nhnt n even ,. , 

= \ (42) 
I t coth -Kt n odd . 

The special case n = was stated without proof by Mehler in 1881 The general 
case is proven in the Appendix. 

Since gn\i) = 1, the final solution becomes a simple superposition of the discrete 
and continuous contributions 

G„ = Gf/^^^ + G™"* (43) 
where the discrete contribution (for n > 2) is 

n/2 L"/2J 



Gd-(z,r) = (^) E - Pi°,")(z)e-("'/2+i/4-((n+i)/2-fc)^). _ 

(44) 



The continuous contribution is 

n/2 



Gr(-,r)=(^i±ij dt^„(t)Pl"(:|,^/,_,,^,(z)e-("V2+i/4+t^). (45) 



with /in(i) given by (^2|). Thus, we arrive at the final result. The probability for 
drawing a given 2x2 matrix Y is 

^(1^,t) = :^ lGo{z,T) + 2Y,Gr.{z,T)cosne\ (46) 
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Figure 1: Plot of G{z,9) for t = 1. Notice the characteristic log-normal tapering of the ridge as 
a function of z, and the nearly Gaussian distribution in 9 around 6 = 0. 

with F{D) given by eqn. ( p^ ) and G'„(z, r) given by eqns. (p3[-^5|). As noted previously, 
the Gn{z,T) are independent of the sign of n so that V is manifestly real. In fig. |^ 
the function G{z, 6) (the expression in parenthesis) is plotted for r = 1. 

6. Conclusions: We have analytically derived the distribution of an infinite 
product of random 2x2 matrices. In statistical image analysis, it may be used as an 
uncommitted prior for morphing and warping Q , with desirable properties not shared 
by the usual priors based on elastic membranes. The distribution of such matrices 
may be evaluated numerically at a moderate cost in computer time and converges 
reasonably fast because of the strong exponential damping. 
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7. Appendix: The Jacobi functions are related to the hypergeometric functions, 

'n+1 



2 ' 2 

with t real, and obey the orthogonality relation 



n + 1 1 — z 

it. — z \- it;l; — - — 



1 + z 



dzP 



(0 



Sjt-t') 



(47) 



(48) 



In order to find /i„(i) for arbitrary n, it is helpful to consider the asymptotic form of 
these functions by using the standard relation for hypergeometric functions 



F (a, b; c; z) 



(1-z) 



T{b)T{c-a) 



F [ a, c ~ b; a ^ b + 1; 



1 



i^-^)-' lifJ"~^l F{b,c-a;b~a+l 



r(a)r(c-6) 



1- z 

1 



This form allows us to see that 

-,(0,n) 



as z ^ oo. Here, 

Ait) = 



P^Li/2+.A^) - 2|A(t)| z-"/^-i/2cos(0, +nnz) 
T{2it) 



r{n/2 + 1/2 + it)T{-n/2 + 1/2 + it) 



2n/2+l/2-it 



(49) 



(50) 



(51) 



and (j)t is the phase of A{t). Using this asymptotic form, we can perform the integral 
in eqn. (48) by using the variable u = log 2, adding a convergence factor of exp (—fiu), 
and finally taking the limit /i — > 0. The result is simply 

2fi 



ix^ + {t-t'Y 



(52) 



The factor in brackets is a familiar representation of 2TT6{t — t') in the limit /x — > 0. 
Standard relations for the gamma function immediately yield eqn. (36). This confirms 
the results of Mehler for the special case 71 = 0. The extension to n > would 
appear to be new. 
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